Numerical Method for Zero- Temperature Vortex-Line Phase Diagrams 






o 

(N 

o 
o 

Oh' 



Welles A. M. Morgado and Gilson Carneiro * 

Instituto de Fisica 

Universidade Federal do Rio de Janeiro 

C.P. 68528 

21945-970, Rio de Janeiro-RJ 

Brazil 

(February 1, 2008) 

A numerical method to calculate equilibrium vortex-line configurations in bulk anisotropic type- 
II superconductors, at zero temperature, placed in an external magnetic field is introduced and 
applied to two physical problems. The method is designed to search for the minimum of the Gibbs 
free-energy in the London approximation and assumes only that the vortex hues are straight and 
arranged in a periodic lattice. Based on these assumptions a simulated annealing algorithm is 
developed to find the vortex- line-lattice unit cell shape, and vortex- line positions within it. This 
algorithm is made fast and accurate by the use of a rapidly converging series to calculate the 
lattice sums entering the vortex-vortex interaction energy. The method's accuracy is illustrated by 
calculating the magnetic induction versus applied field curve for an isotropic superconductor. The 
method is applied to a superconductor with a square lattice of columnar defects to study selected 
regions of the zero-temperature-phase diagram where vortex-line-lattices commensurate with the 
columnar defect lattice exist. 
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I. INTRODUCTION 

A first step for understanding the behavior of type-II 
superconducting systems is to determine the equilibrium 
vortex configurations that occur at temperature zero, 
or the zero-temperature phase diagram. Although this 
problem has been under study for many yearsEJ, there are 
many systems, such as artificial structures, constrained 
geometries and anisotropic superconductors, for which 
some questions remain open. 

The situation of greater interest is that of a supercon- 
ductor subjected to an external magnetic field H. The 
theoretical problem of calculating the zero-temperature 
phase diagram is then to find the vortex configuration 
that minimizes the Gibbs free-energy for a given H. The 
analytical method to solve this problem determines, in- 
stead, the vortex configurations that minimize the energy 
for a given magnetic induction, B, and then calculates 
H from the derivatives of the energy with respect to BEI. 
Analytical solutions are difficult to obtain and are known 
only for a few situations. Numerical methods to solve 
this problem have been proposed. These attempt to find 
the zero-temperature phase diagram of a large number of 
vortices by minimizing the energy, if B is kept-|Cpnstant, 
or the Gibbs free-energy, if H s kept constantclu. These 
methods are severely limited by finite-size effects. 

In this paper we introduce a numerical method to cal- 
culate configurations of straight vortex lines that min- 
imize the Gibbs free-energy of bulk superconductors in 



the London approximation for a given H. The essential 
ingredients of our methods are two. First, we assume 
that the vortex lines are straight and restrict the search 
for the minimum free-energy configurations to periodic 
vortex-line-lattices (VLL) with a given number of vortex 
lines per unit cell. This restriction greatly reduces the 
number of parameters that we have to vary to find the 
minimum. Second, we use a rapidly converging series to 
calculate the lattice sums entering the vortex-vortex in- 
teraction energy in the London approximationa. These 
ingredients allow us to create a simulated annealing algo- 
rithm, using Monte Carlo techniques, that is fast and, as 
a consequence, highly|^ccurate. A closely related method 
has been used in Rcf.El to calculate the zero-temperature 
phase diagram for films under parallel fields. 

To illustrate the accuracy of the method we apply it 
to the well known problem of a bulk superconductor un- 
der a field H to obtain the B x H curve and io compare 
it with known results at high and low fieldgj. We also 
apply the method to find some equilibrium vortex-line 
configurations in bulk superconductors and in films with 
a square lattice of columnar defects. This system is of 
experimental interest because there are several studies of 
vortices in films with a lattice of columnar defects cre- 
ated by lithographic techniques or by bombardment with 
an electron beamEMa. 

This paper is organized as follows. In Sec. |l| we discuss 
the details of the theoretical model, numerical method 
and simulated annealing algorithm. In Sec. fll we report 
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the results of the above mentioned apphcations. Finally 
in Sec. IV we state our conclusions. 



II. THEORETICAL MODEL AND NUMERICAL 
METHOD 

We consider a bulk uniaxially anisotropic superconduc- 
tor with a lattice of columnar defects (CD). We assume 
that both the vortex lines and the CD are straight and 
parallel to the c-axis, and that the vortex lines form a 
periodic VLL. Our problem is then to determine, for a 
given external field parallel to the c-direction, i/, the 
VLL primitive unit cell, the number of vortex lines and 
their positions within this cell, that give the absolute 
minimum of the Gibbs free energy per volume 
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where Ew is the vortex line-vortex line inter- 
action energy, -Ev-cd is the vortex-CD-lattice in- 
teraction energy, B is the magnetic induction, 
Hci = 47re/0o is the lower critical field, and 
e = (0o/167r^A^^) In (Xab/S.) is the vortex hue self-energy 
(Aq6 = penetration depth for currents parallel to the ab- 
plane, ^ — vortex core radius). 

We first solve the simpler problem of finding the VLL 
with a fixed number of vortex lines per unit cell, riy, that 
minimizes G. To find which one of these gives the abso- 
lute minimum of G we have to compare the free energies 
of the VLL obtained for the same H and different n„, if 
these are distinct. 

Assuming that the VLL unit cell is defined by the vec- 
tors Li and L2, as shown in Fig. |l|, or by the correspond- 
ing reciprocal lattice vectors G, the London approxima- 
tion expression for E^ is 
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where Ac is the unit cell area and x^ (a = l,...,n„) 
denotes the vortex lines positions within the unit cell. 
For the VLL shown in Fig. |l|, _B = ni,(^o/^c and Ac = 
L1L2 sin0. 

The vortex-line CD-lattice interaction energy is as- 
sumed to be the sum of the single vortex-single CD in- 
teraction energies, namely 

N n„ 

E^-cd = -TTT- y^ y^ y^ [/v-cd(Rj + x^ - r^) , (s) 



J = l Q=l R-d 



where TV is the number of vortex-lattice unit cells, Rj 
denotes the cell positions, R^ denotes the CD positions 



and Uy-cd is the single vortex- single CD interaction en- 
ergy per unit length. The validity of Eq. (|^) requires 
that the CD radius R be such that R ^ 2^, so that only 
a single vortex line can be trapped by the CD, and th| 
the vortex lines mean separation is large compared to 
We use for t/v-cd an approxiiB|ation based on the London 
theory results derived in Refs, namely 
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This expression interpolates smoothly between the large 
r behavior and the value of the energy per unit length 
of a vortex line pinned by a CD ~ eh\{^/2£^/ R), for 
R^^/2^. 

In order to study both bulk samples and films of thick- 
ness D < \ab in the above described framework, we use 
in the expression for E^^, Eq. (ph, instead of Aah an ef- 
fective penetration depth A ~ ^bl^- This mimics the 
dominant logarithmic vortex-vortex interactions at dis- 
tances short compared to the film screening length. The 
vortex-CD interaction is assumed to be the same for bulk 
and films. 

To numerically minimize G, Eq. (|l|), it is necessary to 
evaluate i?vv, Eq. (0). Because the lattice sum in this 
expression converges very slowly, a numerical method 
based on i1— requires considerable amount of coinputer 
time. DoriaQ showed that the lattice sum in Eq. (g) can 
be transformed into a series that converges much faster. 
We use this series in our calculation to generate an efh- 
cient numerical algorithm. 

Our minimization procedure, based on standard simu- 
lated annealing and Monte Carlo techniquesE2l, is as fol- 
lows. For a given H and n.u, we start with a chosen unit 
cell shape and vortex lines located within it. By moving 
the vortex lines, one at a time, and deforming the unit 
cell continuously we generate, using a Metropolis algo- 
rithm, the equilibrium configurations of the cell shapes 
and vortex lines positions within it, for a given fictitious 
temperature T. The changes in G caused by vortex line 
motion or by unit cell deformation are calculated using 
Doria's fast convergent series for E'w, and a direct sum- 
mation over a large CD-lattice for E^-cd- The VLL that 
minimizes G is identified with the configuration at very 
low T. In order to reduce trapping in metastable states, 
T is cycled appropriately. 

The VLL found by the above described method can be 
either commensurate or incommensurate with the CD- 
lattice. Only the commensurate ones correspond to true 
minima of the Gibbs free-energy, since we restrict our 
search for the G minima to periodic vortex-line arrange- 
ments. The incommensurate VLL, being periodic, feel 
only the space average xif the vortex-CD-lattice poten- 
tial, that is, a constantlld. Consequently, these are tri- 
angular VLL that minimize G, Eq. ™), in the absence 
of the vortex-CD- lattice interaction. As discussed in 



detail in Sec. [II B, we expect that in situations where 



vortex-vortex interactions are strong, such as the ones 
considered here, our method is capable of predicting not 
only the commensurate phases, but also of giving reason- 
able estimates for the range of H values over which they 
minimize G. 



III. APPLICATIONS 

In this section we apply the method described above 
to two physical problems. 

A. Defect-free superconductor 

To illustrate the accuracy of our method, we first apply 
it to the well known problem of vortex lines in a defect- 
free superconductor subjected to an external magnetic 
field H pointing in the c-direction. In this case G is min- 
imized by a triangular VLL with a single vortex line in 
the primitive unit cell. We minimize G for several H as- 
suming Uy = 1 and n^, = 2. We find the triangular lattice 
in all cases considered. For a given H , we compare the 
minima obtained for ny = 1 and ni, = 2. We find for 
riy — 2 a, unit cell that has twice the area as that for 
riv = 1 and a value of G that differs from that for n^ = 1 
by less than 3 parts in 10^. We also find that B, cal- 
culated by our method, is very accurate when compared 
with the known aijialytic expressions for H ^ Hd and for 
ffci < i? <C Hc'^. This is illustrated in Fig. | for the 
M X H curve, instead of the B x H because B ^ H. The 
errors seen in this figure indicate that our B values differ 
from the theoretical ones by less than 1%. 



B. Regular array of columnar defects 

Now we consider a superconductor with a CD-lattice. 

We assume that the CD radius is i? = 2^ = 0.02Aab 
and that the CD-lattice is square with lattice constant 
ttcd = 0.6Aa6. For films we assume that A = 5Xab- 

First we look for the matching phase for which B^ = 
'Po/^cd ^^'^ ^^^ VLL is square and commensurate with 
the CD-lattice, with one vortex-line pinned at each CD. 
We find that for both bulk and films with n,„ = 1 the 
matching phase minimizes G in the range of H values 
H< < H < H>. Wc obtain for bulk H< ^ I.OIB^, 
H> = 1.10^0, and for films H< = H> = l.OOB^. When 
H is just outside this range we find that a VLL incom- 
mensurate with the CD-lattice minimizes G. The H- 
range for films is smaller than that for bulk due to the 
value of A = 5Xab- 

We recall that by our method i/5 and H^ are the 
boundaries in the BxH phase diagram where the match- 
ing state becomes unstable with respect to an incommen- 
surate triangular VLL. We may ask what is the relation- 



ship between these boundaries and true ones obtained if 
G were minimized without the restrictions imposed by 
our method. To answer this question it is necessary to 
guess the state to which the matching phase would be- 
come unstable to at these boundaries. One possibility is 
that the matching state becomes unstable with respect 
to the state obtained from it by the addition of vacan- 
cies or interstitialalj. Another possibility is that at these 
boundaries a cornxaensurate-incommensurate (CI) tran- 
sition takes placeO. We believe that the former is pos- 
sible only if the vortex-CD interaction is stronger than 
vortex-vortex interactions, which is not the case here. In 
our calculation the parameters are such that these two 
interactions are of comparable strength so, we believe, 
that the true boundaries here correspond to a CI transi- 
tion. It is known that for strong vortex-vortex interac- 
tions, and not too close to the CI transition point, the 
incommensurate state differs from the triangular lattice 
by small incommensurate elastic distortionaij. Thus, we 
believe that our method makes a reasonable approxima- 
tion for the incommensurate state (except close to the 
CI transition), and that the estimates obtained from it 
for the matching state phase boundaries are close to the 
true ones. 

We also apply our method with n„ = 2 a value of H 
in the range iJ5 < H < H^ . We find that the matching 
state with a rectangular unit cell twice as large as that 
for 71^ = 1 minimizes G and that the G- value differs from 
that for n„ = 1 by less than 1 part in 10^. 

Next we study VLL commensurate with the CD-lattice 
for B < i?0. These can be characterized as follows. We 
assume that the CD-lattice primitive unit cell is a square 
of side ttcd with unit vectors e^ and ey oriented along 
the sides. We expect that for these VLL all vortex lines 
are pinned by the CD. Thus, the VLL primitive unit cell 
vectors, ai and sl2 can be written as 
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If there are riy vortices in the unit cell, the magnetic 
induction in this state is 



Br — Brf, 



I nij,n2y - niyn2x \ 



(7) 



where n^ is restricted hy n^ <| nixn2y — niyn2x \- Thus, 
when B / B^ is rational, there are many non-equivalent 
VLL satisfying the above stated conditions. We apply 
our numerical method to determine which one of these 
minimize G, and to estimate the range of H values over 
which this minimum is stable. 

Wc study B = B4,/b, 5^/4, B^/S, 5^/2, 25^/3 with 
riy — \ and B — 2B^/5 with «„ — 2. We find that for 
films the commensurate VLL shown in Fig. ^ minimize 
G with H = B. The H range where these minima are 
stable is around if = B for which these VLL minimize 
G is less than 10~^B^. For bulk we find that for Bc/^/A 
and B^/2 the same VLL minimize G for H = 0.3305^ 



and H = Q.528B^, respectively. The iJ-range was not 
investigated in detail in this case. 

Our results for B — B^/A and B ~ B^/2 agree wr 
recent Lorentz microscopy experiments reported in Ref. 
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IV. CONCLUSIONS 

In conclusion then we develop a numerical method to 
calculate the zero-temperature phase diagram of a bulk 
superconductor in the presence of an external field H and 
show that it gives accurate results for the B x H curve 
of an ideal superconductor and that it can make predic- 
tions about the zero-temperature phase diagram of the 
non-trivial problem of vortex lines in the presence of a 
CD-lattice. The method can be extended to study the 
zero-temperature phase diagram for anisotropic super- 
conductors in tilted fields and to long cylinders of rect- 
angular cross sections in axial fields. Work along these 
lines is under way and will be reported elsewhere. 
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FIG. 1. Vortex-line-lattice unit cell with n„ = 3. 





FIG. 2. Magnetization for a defect free superconductor ver- 
sus H compared to theoretical results for low H ( full curve) 
and intermediate H ( dashed curve). 



FIG. 3. Commensurate vortex-line-lattices that minimize 
the Gibbs free-energy. 



